gusucode.com > 阵列信号处理书的源码 > MATALB 程序/13.信源数估计MDL算法MATLAB程序/mdl_sourcenumber.m
clear all, close all, K = 6; % Number of antenna snr = -2; % SNR theta =[10,16,20] % DOA Sample = [10 20 30 40 60 80 100 120 150 200 300 400 600 900 1200]'; % Number of snapshots C = 2; Ntrial = 200; jj = sqrt(-1); % ----------------------------------------------------------------------- % Ndoa = size(theta,2); Nsample = length(Sample); pdf_MDL = zeros(Nsample,Ndoa+2); proposedMethod = zeros(Nsample,Ndoa+2); Num_ref = [0:Ndoa+1]'; for nNsample =1:Nsample number_dEVD = zeros(Ndoa+2,1); SNR = ones(size(theta',1),1) * snr; T = Sample(nNsample); for nTrial = 1:Ntrial %================================================================== source_power = (10.^(SNR./10)); source_amplitude = sqrt(source_power)*ones(1,T); % Source standard deviation. source_wave = sqrt(0.5)*(randn(T,Ndoa) + jj*randn(T,Ndoa)); st = source_amplitude.*source_wave.'; d0 = st(1,:).'; nt = sqrt(0.5)*(randn(K,T)+jj*randn(K,T)); A = exp(jj*pi*[0:K-1]'*sin(theta)); xt = A*st + nt; % received signal %======================= MDL method ============================== [Ke,N]=size(xt); Rx = (xt*xt')./T; [u,s,v] = svd(Rx); sd = diag(s); a = zeros(1,K); for m = 0:K-1 negv = sd(m+1:K); Tsph = mean(negv)/((prod(negv))^(1/(Ke-m))); a(m+1) = T*(K-m)*log(Tsph) + m*(2*K-m)*log(T)/2; end [y,b] = min(a); dEVD = b - 1; p_dEVD = find(dEVD - Num_ref == 0); number_dEVD(p_dEVD) = number_dEVD(p_dEVD) + 1; %---------------------------------------------- end %for nTrial pdf_MDL(nNsample,1:end) = number_dEVD'/nTrial; end % %============================================ figure, semilogx(Sample,pdf_MDL(:,Ndoa+1),'b:*') legend('MDL') ylabel('Probability of Detection') xlabel('Number of Snapshots') axis([Sample(1),Sample(length(Sample)),0,1])